From 5078e8d512c8e695f218595f7939065b02cba955 Mon Sep 17 00:00:00 2001
From: =?utf8?q?=C3=98yvind=20Kol=C3=A5s?= <pippin@gimp.org>
Date: Wed, 16 Nov 2016 17:42:33 +0100
Subject: [PATCH] add single precision versions of linear/gamma conversions

---
 babl/base/pow-24.c | 54 ++++++++++++++++++++++++++++++++++++++++++++++
 babl/base/pow-24.h |  2 ++
 babl/base/util.h   | 15 +++++++++++++
 3 files changed, 71 insertions(+)

diff --git a/babl/base/pow-24.c b/babl/base/pow-24.c
index a3bb36e..03fca43 100644
--- a/babl/base/pow-24.c
+++ b/babl/base/pow-24.c
@@ -71,3 +71,57 @@ babl_pow_1_24 (double x)
     y = (7./6.) * y - z * ((y*y)*(y*y)*(y*y*y));
   return x*y;
 }
+
+//////////////////////////////////////////////
+/* a^b = exp(b*log(a))
+ *
+ * Extracting the exponent from a float gives us an approximate log.
+ * Or better yet, reinterpret the bitpattern of the whole float as an int.
+ *
+ * However, the output values of 12throot vary by less than a factor of 2
+ * over the domain we care about, so we only get log() that way, not exp().
+ *
+ * Approximate exp() with a low-degree polynomial; not exactly equal to the
+ * Taylor series since we're minimizing maximum error over a certain finite
+ * domain. It's not worthwhile to use lots of terms, since Newton's method
+ * has a better convergence rate once you get reasonably close to the answer.
+ */
+static inline float
+init_newtonf (float x, float exponent, float c0, float c1, float c2)
+{
+    int iexp;
+    float y = frexpf(x, &iexp);
+    y = 2*y+(iexp-2);
+    c1 *= M_LN2*exponent;
+    c2 *= M_LN2*M_LN2*exponent*exponent;
+    return y = c0 + c1*y + c2*y*y;
+}
+
+/* Returns x^2.4 == (x*(x^(-1/5)))^3, using Newton's method for x^(-1/5).
+ */
+float
+babl_pow_24f (float x)
+{
+  float y = init_newtonf (x, -1.f/5, 0.9953189663f, 0.9594345146f, 0.6742970332f);
+  int i;
+  for (i = 0; i < 3; i++)
+    y = (1.f+1.f/5)*y - ((1./5)*x*(y*y))*((y*y)*(y*y));
+  x *= y;
+  return x*x*x;
+}
+
+/* Returns x^(1/2.4) == x*((x^(-1/6))^(1/2))^7, using Newton's method for x^(-1/6).
+ */
+float
+babl_pow_1_24f (float x)
+{
+  float y = init_newtonf (x, -1.f/12, 0.9976800269f, 0.9885126933f, 0.5908575383f);
+  int i;
+  float z;
+  x = sqrtf (x);
+  /* newton's method for x^(-1/6) */
+  z = (1.f/6.f) * x;
+  for (i = 0; i < 3; i++)
+    y = (7.f/6.f) * y - z * ((y*y)*(y*y)*(y*y*y));
+  return x*y;
+}
diff --git a/babl/base/pow-24.h b/babl/base/pow-24.h
index e822297..2fb338f 100644
--- a/babl/base/pow-24.h
+++ b/babl/base/pow-24.h
@@ -21,5 +21,7 @@
 
 extern double babl_pow_1_24 (double x);
 extern double babl_pow_24 (double x);
+extern float  babl_pow_1_24f (float x);
+extern float  babl_pow_24f (float x);
 
 #endif
diff --git a/babl/base/util.h b/babl/base/util.h
index cfb17b0..3472523 100644
--- a/babl/base/util.h
+++ b/babl/base/util.h
@@ -78,6 +78,14 @@ babl_linear_to_gamma_2_2 (double value)
     return 1.055 * babl_pow_1_24 (value) - 0.055;
   return 12.92 * value;
 }
+static inline float
+babl_linear_to_gamma_2_2f (float value)
+{
+  if (value > 0.003130804954f)
+    return 1.055f * babl_pow_1_24f (value) - 0.055f;
+  return 12.92f * value;
+}
+
 
 static inline double
 babl_gamma_2_2_to_linear (double value)
@@ -86,6 +94,13 @@ babl_gamma_2_2_to_linear (double value)
     return babl_pow_24 ((value + 0.055) / 1.055);
   return value / 12.92;
 }
+static inline float
+babl_gamma_2_2_to_linearf (float value)
+{
+  if (value > 0.04045f)
+    return babl_pow_24f ((value + 0.055f) / 1.055f);
+  return value / 12.92f;
+}
 
 #else
   #define linear_to_gamma_2_2(value) (pow((value), (1.0F/2.2F)))
-- 
2.30.2